Electron localization in the insulating state: 
Application to crystalline semiconductors 
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We measure electron localization in different materials by means of a "localization tensor", based 
on Berry phases and related quantities. We analyze its properties, and we actually compute such 
tensor from first principles for several tetrahedrally coordinated semiconductors. We discuss the 
trends in our calculated quantity, and we relate our findings to recent work by other authors. We 
also address the "hermaphrodite orbitals", which are localized (Wannier-like) in a given direction, 
and delocalized (Bloch-like) in the two orthogonal directions: our tensor is related to the optimal 
localization of these orbitals. We also prove numerically that the decay of the optimally localized 
hermaphrodite orbitals is exponential. 
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I. INTRODUCTION 

A nonmetal is distinguished from a metal by its vanish- 
ing conductivity at low temperature and low frequency: 
we use here the term "insulator" to include any non- 
metal, like the semiconducting materials which are the 
case studies actually addressed in this work. 

Within classical physics, the qualitative difference be- 
tween an insulator and a metal is attributed to the nature 
of the electronic charge: either "bound" (Lorentz model 
for insulators) or "free" (Drude model for metals). In 
other words, electrons are localized in insulators and delo- 
calized in metals. In a milestone paper published in 1964, 
W. Kohn characterized the insulating state of matter in 
a way which is reminiscent of the classical picture: he 
gave evidence that the main feature determining the in- 
sulating behavior of matter is electron localization in the 
ground-state wavefunction.Q Although this work mainly 
addressed correlated many-electron systems, its message 
is very relevant even for materials where an independent - 
electron description is quite adequate, as the semicon- 
ductor crystals studied here. Recently a novel measure 
of electron localization — different from Kohn's one — was 
proposed by Resta and Sorella,el hereafter cited as RS. 
Their approach-ds, deeply rooted into the modern theory 
of polarizationflEI 

Metals and insulators reveal their qualitative difference 
when static dielectric polarization is addressed. Sup- 
pose we expose a finite macroscopic sample to an electric 
field, say inserting it in a charged capacitor. In metals 
polarization is trivial: universal, material-independent, 
due to surface phenomena only (screening by free car- 
riers). Therefore polarization in metals is not a bulk 
phenomenon. The opposite is true for insulators: macro- 
scopic polarization is a nontrivial, material-dependent, 
bulk phenomenon. We can therefore phenomenologically 
characterize an insulator, in very general terms, as a 
material whose ground wavefunction sustains a nonzero 
bulk macroscopic polarization whenever the electronic 
Hamiltonian is non centrosymmctric. If the Hamilto- 



nian is instead centrosymmetric, the polarization van- 
ishes but remains a well defined bulk property, at vari- 
ance with the metallic case. The phenomenological link 
between macroscopic polarization and insulating behav- 
ior was first pointed out and exploited — taking advan- 
tage of the modern theory of polarizationBEJ — by RS in 
1999. This approach is based on Berry phases and re- 
lated concepts.Ell Even the RS paper, like Kohn's 1964 
one, mostly concerns correlated systems. Furthermore, in 
order to keep the presentation simple and concise, most 
results are explicitly shown in one dimension, while the 
d-dimensional formulation is only sketched in the final 
paragraphs of RS. In the present paper we provide more 
details on how the RS theory of localization works in 
three dimensions, specializing to a system of nonintcr- 
acting electrons, like the band insulators chosen as case 
studies here. 

Some other important papers must be mentioned at 
this point. In 1997 Marzari and Vanderbilt JlJ here- 
after cited as MV, while not addressing metals at all 
(and hence their difference from insulators), establish 
nonetheless some results which are relevant to the present 
viewpoint. In a very recent comprehensive papeitj 
Souza, Wilkens, and Martin — hereafter cited as SWM — 
generalize and extend in various ways the main finding of 
RS: we adopt here some of their notations. Finally, after 
this work was completed, we became aware of Ref. [l^, 
whose conclusions bear some implications for our results 
shown in Sect. VI. 

The paper is organized as follows. In Sect. || we 
define the basic ingredients providing both polarization 
and localization, namely the expectation values of the 

many-body phase operators zjy for the three Carte- 
sian coordinates, Eq. fl). In Sect. Ill, following RS, we 



show how the modulus of zffl defines a very meaning- 
ful quantity, the localization tensor, for which we adopt 
the SWM notations: such tensor is finite in insulators 



and diverges in metals. In Sect. IV we discuss the main 
properties of the localization tensor, and in Sect. [v| we 
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present first-principle calculations for several elemental 
and binary semiconductors: the main trends are ana- 
lyzed. In Sect. [y] we calculate orbitals which are opti- 
mally localized in a given direction, and whose average 
quadratic spread coincides with the localization tensor. 
We also heuristically check the exponential localization 
of these orbit als, which we call "hermaphrodite orbitals". 
In Sect. VII wc draw our main conclusions. In the Ap- 



pendix we consider a molecule or a cluster an we discuss 
our localization tensor therein, showing ite-.relationship 
to some results of Boys4pcalization theory,Ej well known 
in quantum chemistry.E-51 



II. MANY-BODY PHASE OPERATORS 

We are addressing here, as it is done by MV, a crys- 
talline system of independent-electrons, having in mind a 
Kohn-Sham scheme. The properties of interest, namely, 
macroscopic polarization and electron localization, are 
not properties of the individual KS orbitals: instead, 
they are global properties of the occupied KS manifold. 
As shown in Refs. ^jj^, it proves formally convenient to 
deal with a many-body wavefunction obtained as a 
Slater determinant of occupied orbitals. This determi- 
nant is uniquely determined by the manifold of the oc- 
cupied orbitals and is invariant by unitary transforma- 
tion of these orbitals among themselves: for instance, in 
insulating crystals, an important transformation of this 
class converts the occupied Bloch orbitals into Wannier 
functions J13 Quantities which can be expressed solely in 
terms of are invariant in form under such transforma- 
tions. 

Throughout this 
work — with the exception of the Appendix — we adopt 
periodic Born-von-Karman boundary conditions (BvK) 
on a large cell, multiple of the crystalline elementary cell. 
The quantities of interest are intensive and have a well 
defined thermodynamic limit, while the wavefunction it- 
self becomes an ill-defined mathematical object in that 
limit. 

For the sake of simplicity, we assume a simple cubic 
cell of side a and a large BvK cell of side L — Ma. 
More general structures can be dealt with using scaling, 
similarly to what shown e.g. in Ref. ^ or in SWM. 
The thermodynamic limit corresponds to M — > oo, 
while practical calculations are performed at finite, and 
possibly large, M values. The spinorbitals ip (spin-up) 
and -0 (spin-down) may be chosen of the Bloch form. In 
the finite system there are M 3 allowed Bloch vectors q s , 
arranged on a regular mesh in the unit reciprocal cell, 
where s = (si, S2, S3) and 

27T 

q. s = jj-{si,s 2 ,s 3 ), s a = 0,l,...,M-l. (1) 

We adopt a plane-wave-like normalization for the Bloch 
orbitals: 



■73 / dr Cq s ( r )VVq s , (r) = 8 nn .8 

L JBvK cell 



ss' ; 



(2) 



If the system is insulating with rib doubly occupied 
bands, there are iV = 2ni,M 3 independent spinorbitals, 
out of which we write a single-determinant many-body 
wavefunction for N electrons: 



(3) 



where the product runs over all occupied bands and all 
mesh points, A is the antisymmetrizer operator, and 
the factor ensures that the TV-body wavefunction is 
normalized to one on the hypercube of side L. If, instead, 
the system is metallic, then the many-body wavefunction 
^ can still be written in the form of Eq. (||), but where 
not all the Bloch vectors of a given band are included in 
the product. 

According to Refs. [|,|[ the key quantities to deal with 
both macroscopic polarization and electron localization 
are expectation values of "many-body phase operators" . 
For a three-dimensional system there are three such 
operators, one for each Cartesian direction. We indicate 
as , where a is a Cartesian index, their ground-state 
expectation values: 



(x) 



(*|C 



(4) 



and analogously for y and z directions. This remarkably 
compact expression is very general and applies as it 
stands even to correlated and/or disordered systems: 
here we specialize to a crystalline system of independent 
electrons, whose wavefunction ^ assumes the form of 
Eq. (||), where the product indices have to be differently 
specified in the insulating and metallic cases. 



We may conveniently recast zff as an overlap 



(x) 
Z N 



<*!*>, 



(5) 



where \& is the Slater determinant of a different set of 
Bloch spinorbitals: 



0nq 3 (r) = e l M°- X 4> n q s ( r ), 



(6) 



and analogously for the bar (spin-down) ones. According 
to a well known theorem, the overlap between two single 
determinant wavefunctions is equal to the determinant 
of the N x N overlap matrix built out of the occupied 
spinorbitals. Since the overlaps between different-spin 
spinorbitals vanish, and those between equal-spin ones 
are identical in pairs, we can write 



Z N 



(det Sf 



(7) 



where S is the overlap matrix between spatial orbitals, 
having size N/2 — n^A/ 3 . Its elements are: 
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S„ qs ,„< qs , = 73 / dr V4j s ( r )^n'q s ,(r) 

JBvK cell 

= 4/ rfr< qs (r)v q „(r)e t ^ I+ ^'"- t ""-», (8) 

-° JBvK cell 

where the u's are the periodic functions in the Bloch 
orbitals. 

The matrix S is very sparse: in fact, given the 
geometry of the q s 's on the regular reciprocal mesh (see 
Eq. (0) , the overlap integrals in Eq. (||) are nonvanishing 
only if si = s[ + 1, S2 = s' 2 , and S3 — s' 3 . We express the 
nonvanishing elements in terms of a small overlap matrix 
S, of size rib X i&: 

Snn'i^q') = (w nq |w„/ q ,) = i I dru* q (r)tt n / q /(r). (9) 

a J cell 

Owing to the sparseness of S, its determinant factors into 
products of determinants of small matrices S. 

In the insulating case we use the wavefunction of 
Eq. (H), where all the Bloch vectors of a given band are 
occupied: the factorization is then 

z^ )1/2 = det S = TJdet S(q sl+ i, S2tSa , q Sl , S2>S3 ). (10) 

S 

We get a more compact notation upon defining 

Aq(*> = f (1,0,0) = ^-(1,0,0), (11) 

which is the vector connecting nearest neighbor q points 
in the x direction. We have then: 

z^^UdctS^ + A^,^); (12) 

s 

In the metallic case, instead, the zj^'s are identically 
zero. This is easily understood by looking at the simple 
case with only one band. Suppressing the band index 
the small overlap matrix becomes a c-number S^q, q'), 
and Eq. ( JT~0| ) becomes a product of c-numbers, with no 
determinant to evaluate. In an insulator this product 
runs over the whole q s mesh, and all factors are nonzero; 
in a metal the analogous product runs only on the q s 's 
within the Fermi surface. Looking at the definition of 
^(QjQOj Eqs. (g) and (^|), it is clear that there exists 
at least one occupied q s , adjacent to the Fermi surface, 
such that jS^q^qs') vanishes for all occupied s' . This is 
enough to imply that vanishes as well. 

The complex numbers zff are ground-state expecta- 
tion values, and do not access any spectral information. 
Yet they qualitatively discriminate between insulators 
and metals: they are in fact nonvanishing in the former 
materials, and vanishing in the latter ones. This shows, 
according to RS, that there is a qualitative difference in 
the organization of the electrons in their ground state. 
It is remarkable that, in the present case, such difference 
shows up already at finite N, before the thermodynamic 
limit is taken. 



III. ELECTRON LOCALIZATION 

In centrosymmetric materials the expectation values 
Zjq are real (provided the origin is chosen at a cen- 
trosymmetric site), while in noncentrosymmetric mate- 
rials they are in general complex: their phases define 
then the Cartesian components of the macroscopic po- 
larization in suitable units.Qtl In the metallic case the 
2jv s vanish and the polarization is ill defined, in agree- 
ment with the phenomenological viewpoint illustrated in 
Sect. Q. We address electron localization using the mod- 
uli of these same z^s. Following RS, electron local- 
ization is measured by a squared localization length in 
one dimension, and by a "localization tensor" in three 
dimensions. This tensor is an intensive quantity, has the 
dimensions of a squared length, and measures the local- 
ization of the many-electron system as a whole: in the 
present case, it is a global property of the occupied KS 
manifold. The localization tensor is finite for insulators 
and diverges for metals. _ 

In the very recent SWM paperlij it is shown, among 
other things, that the RS localization tensor is related 
to the mean-square quantum fluctuation of the polariza- 
tion: it is a second cumulant moment, which can be very 
elegantly extracted from a moment generating function. 
We adopt throughout notations inspired by SWM, and 
we indicate the localization tensor as (r a rp) c , where the 
subscript stays for "cumulant" . For a material having 
cubic or tetrahedral symmetry, like the semiconductors 
considered in the present case studies, the localization 
tensor is isotropic: its only independent element is {x 2 ) c . 
Its expression is provided by RS, whose Eq. (18) we re- 
cast here as 

= *\*&\*> (13) 

and the thermodynamic limit is understood. For a 

(x) 

metal z N vanishes and the localization tensor is formally 
infinite, even at finite N. For an insulator, whose 
wavefunction has the form of Eq. (||), we get from 
Eq. © 

\zP\ = Y[det 5t(q 5 ,q s +AqW)5(q s ,q s +AqW); (14) 

S 

«< = -(£) 3 2^»! 2 - < 15 > 

Eqs. ( [T4] ) and (|l^) are the typical expressions im- 
plemented in our test-case calculations discussed be- 
low. The thermodynamic limit is obtained as usual for 
M — > 00 and takes, not surprisingly, the form of an in- 
tegral performed over the reciprocal unit cell, or equiva- 
lently over the first Brillouin zone. The integral is: 




3 



'dq x 



l Jiq/ 



(16) 



The proof is relatively straightforward, starting from 
Eq. (|l6|) and discretizing integrals and derivatives on the 
mesh defined in Eq. (jl]). 

Expressions such as Eq. (|l^) and similar ones had 
appeared in the literature before,t£l in relationship to 
Wannicr functions. By means of an expression of this 
kind, MV define a ground-state quantity Vl\ which sets 
a lower bound for the second (spherical) moments of the 
Wannier functionsJiil More precisely, for an insulator with 
rib occupied bands (hence rib Wannier functions per cell) 
such second moment is no smaller in average than fli/ rib- 
It is worth mentioning at this point that the logic of the 
MV paper goes backwards with respect to the present 
approach: first they provide a continuum theory, and 
then they discretize for computational purposes. Their 
discretization is different from Eq. (15), which emerges 
naturally from the present formulation starting from 
the remarkably compact Eq. (|l3|). Both discretizations 
obviously converge to the same M — > oo limit: their 
convergence properties are different, though. 

Specializing MV to a cubic material, RS have found 
the simple relationship fl\ = 3nb{x 2 } c : notice that (ir 2 ) c 
is intensive, while f2i is not such. Building upon MV's 
work, we are now ready to generalize the localization 
tensor to materials of arbitrary symmetry: 



{r a ri3)c 



n fc (27r) : 



d , d 



^(Unq_\-g^-U n/q )( — U nlrl \u n q) J , (17) 



where V c is the cell volume. Notice that the imaginary 
part of the integrand in Eq. (p"7j), being antisymmetric 
in q, cancels in the integral, such that the localization 
tensor is real. Even the offdiagonal elements, as defined 
in Eq. (|l7|), have a finite-iV counterpart in terms of 
many-body phase operators. 

For an insulating crystal of arbitrary symmetry, fli 
as defined by MV equals rib times the trace of our 
localization tensor {r a rp) c . In a metal, expressions like 
Eqs. (|l6|) and ([l7]) do not make much sense, consistently 
with the fact that our finite- N expression, Eq. (13), is 
formally infinite at any TV value. 



IV. PROPERTIES OF THE LOCALIZATION 
TENSOR 

We have already emphasized that the localization 
tensor is a property of the occupied KS manifold as a 
whole. The main quantity which defines such manifold 
is the (spin-integrated) single-particle density matrix 



p, which coincides with twice the projector P over the 
occupied KS orbitals: this projector is invariant by 
unitary transformations of the orbitals. Using Bloch 
eigenf unctions the projector reads, for an insulator with 
rib occupied bands: 

P(r,r') = lp(r,r') = ^3 £ / ^ ^ q (r)C q (r')- 

(18) 

The localization tensor (in the thermodynamic limit) has 
been written as a Brillouin-zone integral in Eq. (|l7|). 
This integral can be identically transformed into a par- 
ticularly simple expression whose only ingredient is P: 

{r a rp) c = ^- f dr f dr' (r - r%(r - v') [3 |P(r, r')| 2 , 

^ n b Jccll Jail space 

(19) 

which is the second moment of the (squared) density 
matrix in the coordinate r — r'. The proof of the 
equivalence between Eq. ( |l9| ) and Eq. (17) can be worked 
out using the same algebra appearing in Ref. |l6|: for 
a different argument proving the same result, see the 
Appendix. 

We have arrived at Eq. (|l^) considering an insulat- 
ing crystali-sOjiar. In this case we know, under general 
argumentsJilrE] that P(r, r') is asymptotically exponen- 
tial in the argument |r — r'|: this confirms that the in- 
tegral over all space in Eq. (|l^) converges and the lo- 
calization tensor is therefore finite. At this point, it is 
worthwhile to apply the general form of Eq. (|l^) to the 
metallic case. For the simplest metal of all, the-iree elec- 
tron gas, the density matrix is known exactly £il 



P(r,r') = -p(r,r') 



3n ji(fc F |r - r'|) 
2 fc F |r-r'| ' 



(20) 



Replacement of Eq. (|20|) into Eq. dl9|) results in a diverg- 
ing integral, thus confirming that our localization tensor 
is formally infinite in this paradigmatic metal. Other, 
more realistic, metals feature this same divergence. 

The fact that the density matrix p(r, r') is short- 
range in the variable p-,— r' has been named "nearsight- 
edness" by W. Kohn.Ej The second moment expression 
in Eq. (fL9|) shows that our localization tensor is indeed 
a meaningful quantitative measure of such nearsighted- 
ness. We are going to analyze below the major trends 
over an important class of materials: tetrahedral semi- 
conductors. We mention at this point that a conceptu- 
ally different measure of the nearsightedness of a given 
electronic ground state focuses instead on the exponent 
governing the exponential decay of p(r, r') in insulators: 
some case studies have been recently investigated.E3o 

We have already observed that some of our findings 
are closely related to the previous work by MV. These 
authors' main interest were the "optimally localized 
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Wannier functions", i.e. those localized orbitals which 
minimize the average spherical moment. They prove, 
among other things, that such moment is strictly larger 
than the trace of our localization tensor. Building on 
their results, it is straightforward to attribute a similar 
meaning to the tensor itself: for any transformation of 
the occupied orbitals into a set of unitarily equivalent 
ones, the second moment in a given direction can be no 
smaller than the localization tensor, projected on that 
direction. 

Since we are going to apply our results to cubic 
materials only, we focus on those orbitals which minimize 
in average the quadratic spread (second moment) in the x 
coordinate. The present formalism makes the definition 
of these orbitals particularly simple: they are in fact the 
eigenfunctions of the position operator x, projected over 
the occupied manifold. Calling 5 = PxP this operator, 
its expression in the Schrodinger representation is: 

S(r,r')= ( dr" P{y,y")x"P{y",y'). (21) 

J all space 

Notice that x is incompatible with BvK boundary con- 
ditions and its matrix elements over Bloch states are ill 
defined; nonetheless, S is — in insulators — a well defined 
operator, which maps any vector of the occupied mani- 
fold into another vector of the same manifold. This fact 
owes to the exponential localization of P in Eq. (^l|). 

The relationship between 5 and the orbitals optimally 
localized in the x direction is easily proved borrowing 
some results from MV; for a different argument leading to 
the same proof, see the Appendix. We also notice an im- 
portant difference with respect to the three-dimensional 
localization explicitly considered by MV. While the trace 
of the localization tensor provides a lower bound for 
three-dimensional localization, its element (x 2 ) c provides 
instead a genuine minimum for one-dimensional local- 
ization (in a cubic material). This qualitative difference 
owes to the fact that, while one can manifestly diagonal- 
ize PxP, one cannot simultaneously diagonalize PxP, 
PyP, and PzP. 

We end this Section about the properties of the 
localization tensor with a most important issue: is 
(r a r/3)c a measurable quantity? The answer, due to 
SWM, is "yes" . They prove the identity: 



(r Q r /3 ) c 



hV c 
2ne 2 nb 



dio 



Re a a p(uj), 



(22) 



where a a p is the conductivity tensor. Notice that the left 
hand side, as emphasized throughout the present work, is 
a property of the electronic ground state, while the right 
hand side is a measurable property related to electronic 
excitations: therefore Eq. (|22j) must be regarded as a 
sum rule. The frequency integral in Eq. (j2^) diverges in 
metals and is finite in insulators, as obviously expected. 
Since in the latter materials there is a gap for electronic 
excitations, Eq. (|22|) immediately leads to the inequality: 



{ra r \ < Wc 

13 ° 2TTe 2 n b £ g Ju 



p CO 

/ dui Re o- a fs(uj) 
Jo 



(23) 



where e g is the direct gap. Using then the oscillator- 
strength sum rule, Eq. ( p3| ) for a cubic material is cast 
as: 



(X* c < 



2m e Sg 



(24) 



Below, we investigate the trends in both members of this 
inequality for our test-case materials. 



V. CALCULATED LOCALIZATION TENSORS 

We have studied several tetrahedrally coordinated 
crystalline materials, from the group IV, III-V, and II- 
VI, having the diamond and zincblende structure. The 
first-principle calculations have been performed within 
density-functional theory in the. local-density approxi- 
mation, using pseudopotentialsEj and plane waves. We 
implement a trivial extension of the formulas presented 
above, using a rectangular unit cell instead of a simple 
cubic one: we thus describe the diamond and zincblende 
structures by means of a tetragonal cell with a lattice 
constant a in the basal plane and c = \/2a. There are 
four atoms per unit cell, whose projections on the c axis 
are equispaced; for the sake of consistency with the for- 
mal results, we take x along the c axis and yz in the 
basal plane. We then use a BvK cell of sides M x c, M y a, 
and M z a, corresponding to a mesh in reciprocal space 
with M x , My, M z points: this allows an easier control of 
convergence. 

We start evaluating at the mesh points the Hermitian 
matrices 



A s = fit ( qs , qs + Aq(-) )S(q s , q s + Aq<*> ) ; 



(25) 



then Eqs. ([14]) and ( |l5| ) are written as 

c 2 M x 

M y M z ^ ^ \ ~4n 2 n b 

" S 2 = l S 3 = l 



M v Mz / 2 M x 

81=1 



(26) 



In general convergence is fast in M y , M z , and slower in 
M x . The expression in parenthesis in Eq. ( p6| ) is precisely 
the one-dimensional expression discussed in detail by RS, 
and the three dimensional one simply obtains from it as 
an average in the (q y , q z ) plane. 

First we show in Fig. |l| the convergence of our expres- 
sions over a genuinely cubic grid, which coincides with 
the one used by MV in their evaluation of the quantity 
fli = 3nb(x 2 ) c . They use a different discretization of 
the same k space integral: both calculations converge to 
the same localization tensor, although our discretization, 
based on Eqs. (14) and (15), converges faster. All the fol- 
lowing results have been obtained using noncubic grids, 
as in Eq. (Ejq), in order to achieve faster convergence. 
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We have systematically calculated well converged lo- 
calization tensors for several elemental and binary semi- 
conductors. In Fig. H we plot the localization tensors ver- 
sus the right-hand member of the inequality in Eq. (24), 
where for the gap e g we have used both (i) the theoretical 
and (ii) the experimental values. In case (i) the inequality 
owes to an exact sum rule and must be satisfied: we are 
therefore checking the internal consistency of the calcula- 
tions. Also, it may be noticed that the inequality is very 
strongly verified. In case (ii) there is no a-priori guaran- 
tee that the inequality is verified, particularly given the 
fact that the experimental gap is systematically larger 
than the KS one. Nonetheless the localization tensor is 
obtained here as a pure ground state property, and it is 
well known that density-functional theory in the local- 
density approximation provides a good representation of 
the ground state, though not of the excitations. a It is 
therefore interesting to verify that even for case (ii) the 
inequality in Eq. (|2j) is strongly verified for all the ma- 
terials considered. 

The localization tensor ranges roughly between 1 and 
3 bohr 2 for all the materials considered, diamond being 
the most localized and germanium the most delocalized. 
The trend is qualitatively expected, in agreement with 
SWM's statement that "the larger the gap, the more 
localized the electrons are". However this is a trend 
more than a strict rule, and indeed a few materials 
show irregularities. Better trends are obtained when 
comparing families of materials: either isoelectronic 
series or isovalent series. In order to enhance such 
regularities, we have heuristically tried a few different 
laws. In Fig. [| we plot our localization tensor versus 
l/sg, using the minimum gaps instead of the direct ones: 
here monotonical trends are very perspicuous. 



VI. MAXIMALLY LOCALIZED 
HERMAPHRODITE ORBITALS 

We actually perform localizing transformations on the 
Bloch orbitals. At variance with the most standard ap- 
proach, we focus on orbitals which are localized in one di- 
rection only, say x, while they arc completely delocalized 
in the yz directions. By analogy with the standard the- 
ory of Wannier functionsjlj one obtains such orbitals by 
integration of the Bloch ones over one component only of 
the Bloch vector. The resulting orbitals are Wannier-like 
in one direction and Bloch-like in the other two: they can 
be therefore called "hermaphrodite orbitals" . Because of 
the same reasons as for ordinary Wannier functions, such 
hermaphrodite orbitals are nonunique: we focus here on 
those hermaphrodite orbitals which are optimally local- 
ized in the x direction. It has been shown above that, 
in the thermodynamic limit, these orbitals are eigenfunc- 
tions of the operator 5, Eq. (|2l|), and their centers are 
the corresponding eigenvalues. It is expedient to consider 
the modified operator 



S(r,r') = / dr" P(r,r")e l " 

J all space 



P(r",r'), (27) 



which to leading order in 1/M X has the same eigenfunc- 
tions as S, and simply related eigenvalues. 

When considering a finite sample with BvK boundary 
conditions — or equivalently a discrete grid in the recip- 
rocal unit cell — the operator S as in Eq. ( pl|) is useless 
because the operator x therein becomes ill defined. In- 
stead the operator EE is well defined, provided the value of 
M x is consistent with the choice of the grid. The integral 
in Eq. ( |27j ) is now performed over the BvK cell and not 
over all space; the projector projects over the finite occu- 
pied manifold, having dimension nbM x M y M z . Choosing 
the Bloch functions as the basis in the occupied manifold, 
the matrix elements of EE are nothing else than the matrix 
S defined in Eq. (||). Therefore in the discrete case the 
hermaphrodite orbitals which achieve optimal localiza- 
tion in the x direction are simply obtained by diagonal- 
izing the matrix S. Since, as already observed, the matrix 
is already diagonal in s 2 and S3, the problem is reduced 
to M y M z independent diagonalizations of submatrices of 
size ribM x . We characterize our orbitals as Wj,s 2 ,s 3 , where 
(*2, S3) is a two-dimensional Bloch-index and j is a one- 
dimensional Wannier-like index, with 1 < j < ribM x . 

We are going to verify that these orbitals indeed 
minimize the average quadratic spread in one-dimension. 
If we define 



dr K S2jS3 (r)|V 



BvK cell 



(28) 



then according to RS the quadratic spread of one given 
hermaphrodite orbital is: 



2 c 2 Ml 

\ W j,S2,S3\ X 1^,82,83/0 = ^ 2 ^ \ Z 3; S 2,S3 I 



(29) 



Taking now the average over all orbitals, and calling 
this quantity \ xx , we get 

I M * M * / C 2 M n » M * 



M y M z ^ ^ \ 4ir 2 n b 

S 2 =l S 3 = l \ ] = 1 



(30) 



We then notice that Wj tS2iS3 are the eigenvectors of S, 
hence the expectation values Zj. S2 . S3 are the correspond- 
ing eigenvalues. Since the product of the eigenvalues 
equals the determinant, standard manipulations prove 
that the average spread A^^squals indeed the lower 
bound (x 2 ) c as given in Eq. (Eq). 

There is a subtlety about the diagonalization of the 
submatrices of S, which are the projection over a cer- 
tain finite dimensional manifold of the operator e l M ^" x . 
Although the operator is unitary, its projection is not a 
unitary matrix, hence the eigenvectors are not exactly 
orthogonal, as instead honest localized orbitals must be: 
this is not a serious problem. In fact the larger is M x , 
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the closer to unitarity the matrix becomes: we know that 
the modulus of its determinant differs by one for a term 
of the order 1 /M x , hence the modulus of each eigenvalue 
differs by one for a term of the order . We recover 

exact orthonormality in the thermodynamic limit; in our 
calculations already at M x ~ 20 deviations from orthog- 
onality are hardly noticeable. 

We have calculated the orbitals Wj. S2 . S3 for several 
crystalline semiconductors: to the purpose of display, 
we call the yz average of \wj, S2 , S3 \ 2 as n\ oc {x), where 
the indices remain implicit. At fixed (s2,Ss) we have, 
given our double cell, 8M X orbitals centered on a BvK 
period of length M x c. There are however at most four 
different shapes, and one obtains all the functions upon 
translations in the x direction (by multiples of c/2) of 
the four basic ones: this is not surprising, since the 
genuine unit cell is one half of our computational one. We 
find that the different shapes are actually always four, 
with the only exception of an elemental semiconductor 
at «2 = S3 = 0. In this very special case the different 
shapes are only two, the orbitals are centered at the 
bond center, and their densities are centrosymmetric: the 
corresponding functions n\ oc {x) are shown in Fig. ^ for 
the case of Si. The most general case is exemplified by 
Fig. [|: it shows the four different n\ oc (x) for the case 
of GaAs, again at s 2 — S3 = 0. None of these four 
w orbitals is therefore centered at a symmetry site, and 
none is centrosymmetrical, although they are obviously 
symmetrically related to each other. About the actual 
value of the quadratic spread (in the x direction) of each 
of the w's, we have found as a general feature that the 
least localized ones are those for S2 = S3 = 0, i.e. at the 
r point in the two-dimensional reciprocal space. 

We now address the long standing issue of exponential 
localization. Exact general results only exist for the gen- 
uine onCpdimcnsional case, where W. Kohn has proved 
long agoEZI that the Wannier functions which minimize 
the quadratic spread (i.e. are optimally localized) have 
an asymptotic exponential behavior. After the present 
work was completed, we became aware of Ref. [l3[ where 
the asymptotic exponential is shown to have a power- 
law prefactor. In three dimensions the problem is utu 
solved, with the exception of some very special cases.EEl 
It has been conjectured by MV that their optimally lo- 
calized Wannier functions enjoy three-dimensional expo- 
nential localization: an analytical proof looks very hard. 
Our hermaphrodite orbitals w are optimally localized in 
the x direction, and therefore in a sense they have one- 
dimensional character: nonetheless, they are have a gen- 
uine dependence on all three coordinates. Therefore an 
analytic proof along the lines of Refs. [l3] and [Tt] is not 
easily extended to our case. Instead, it is simple to use 
the same arguments as given in Ref. 25 in order to prove 
that our w orbitals decay in the x direction faster that 
any inverse power of x. 

Our very elongated BvK cells allow us to study the 
asymptotic behaviour heuristically on our calculated w's. 
The nearest periodic replica of a given w orbital is 



centered at a distance of M x c from it: therefore the 
interesting "asymptotic region" is accessible up to a 
distance somewhat smaller than M x c/2, as it clearly 
appears from Fig. ^. The quantity of choice in order 
to "blow up" the exponential behavior is obviously the 
logarithm of n\ oc (x), which we plot in Fig. ^ (thin solid 
lines) for the case of Ge and for two different (s 2 , S3). It 
is seen that there is a wide region where the plots have 
an overall linear behavior, with superimposed oscillations 
having the crystal periodicity along the x direction (c/2 
in the present case). The slopes at different (s 2 ,s 3 ) are 
very different, though; the n\ oc (x) with the slowest decay 
corresponds to s 2 = S3 = and therefore to the least 
localized, as we previously observed. Next, we filter 
the disturbing periodic oscillations using our favorite 
tool of the macroscopic average.c3 We tried both ways: 
filtering n\ oc (x) first and then taking the logarithm, or 
filtering lnn\ oc (x): the latter turns out to work best. 
The macroscopic filtering is also shown in Fig. ^| (thick 
solid lines): it is easily realized, expecially looking at 
the magnified plot in the lower panel, that there is a 
sizeable region, spanning several cells, where the plotted 
function looks accurately linear with x, hence n\ oc (x) oc 
exp(±6x). We therefore demonstrate "experimentally" 
the exponential localization of our w orbitals. After we 
became aware of Ref. 13, we checked that the power- 
law prefactor suggested therein does not improve the 
quality of our fits. It is hard to assess whether this 
is due to a basic difference between our case and a 
genuine one-dimensional one, or to the limited resolution 
achievable in our selfconsistent three-dimensional finite 
size calculation. Finally, in Fig. we display some 
correlations between the localization length and the 
exponential decay length 1/6 averaged over the two- 
dimensional mesh (s2,S3). 



VII. CONCLUSIONS 

In the present work we provide the three-dimensional 
formulation of the RS theory of electron localization,!!! 
specializing it to the case of independent KS electrons; 
we discjuss-some of its relationships to the MV and SWM 
papers and also (in the Appendix) how-, it relates 
to Boys theory of localization in molecules. E-3 We then 
implement the theory to several materials in the class 
of tetrahedrally coordinated semiconductors. Among the 
results, we find that in general the calculated localization 
length is a monotonical function of the gap, although 
a few materials show irregularities. The trend is more 
regular within a given family (isoelectronic or isovalcnt). 
Finally, we heuristically show that the orbitals which are 
optimally localized in a given direction ( "hermaphrodite 
orbitals" ) show exponential localization. 
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APPENDIX: RELATIONSHIP TO BOYS 
LOCALIZATION IN MOLECULES 

We abandon here the BvK boundary conditions used 
throughout this work, and we consider an iV-electron 
system which is bounded in space. Both the orbitals and 
the many-body wavefunction W are therefore exponen- 
tially vanishing at large distances. Supposing that N is 
even and the state is a singlet, for independent particles 
the wavefunction is the Slater determinant: 



where "Tr" indicates the trace on the electronic coordi- 



1 



Pl<Pl<P2<P2 ■ ■ ■ <PN/2<Pn/2 



(Al) 



The orbitals enjoy no specific symmetry. Any unitary 
transformation of the orbitals produces the same many- 
body ground state (modulo an overall phase): a specific 
choice of the orbitals will be referred to as "choice of 
the gauge" in the following. Obviously all ground state 
properties are gauge-invariant. The density matrix is 
twice the projector over the occupied orbitals: 



P(r, 



2P(r,r') 



N/2 

2£>(i>*(r')- 

i=i 



(A2) 



We are interested in exploiting the gauge freedom in 
order tcv-express the ground state in terms of localized 
orbitals .E£l The standard Boys localization^ in molecules 
consists in minimizing spherical second moments, in 
perfect analogy with MV, which can be regarded as 
the solid-state analogue of Boys localization. Here 
instead we are mostly interested in localizing in one given 
direction, say x. 

For any given choice of the single-particle orbitals <fi, 
the average quadratic spread in the x direction is by 
definition: 

9 N/2 

A L = ^E(^I^ 2 I^)"^N^) 2 )- ( A3 ) 

i=l 

We recast this identically as: 



\2 _ _ 

11 TV 



^2(ipi\x (i - 53 IviKvil ) x \<Pt) 

i 3 

2 



N 



53 K^M'ft)! 2 



(A4) 



The first term in Eq. (A4) is gauge invariant, since we 
can identically write: 

AL = |pTr xPx(l -P) + lj2 \{<Pi\*\<Pj)\ 2 , (A5) 



nate. The gauge-invariant term in Eq. (A5) can be re- 
garded as the xx element of a more general tensor, which 
turns out to be the molecular analogue of our localization 
tensor. We use the same notation for molecules and for 
crystals: 



(r a rp) c = — Tr r a Prp{\ 



(A6) 



If we look for the orbitals which minimize the average 
spread in the x direction, the solution, after Eq. (|A5|), 
is provided by those orbitals which diagonalize the po- 
sition operator x, projected over the occupied manifold. 
Obviously, a set of orthonormal orbitals which diagonal- 
ize PxP can always be found, since PxP is a Hermitian 
operator. The quadratic spread of the se o rbitals equals 



3 , the gauge-invariant part in Eq. (AE). If we are in 



terested instead in minimizing the spherical second mo- 
ment, in general we cannot diagonalize simultaneously 
PxP, PyP, and PzP. Therefore the spherical spread 
will be in general strictly larger than the Cartesian trace 
of the localization tensor. This is a key feature in the 
work of Boys,Efl and MV as well. 

We have defined the localization tensor in Eq. ( [AC] ). 
With an obvious generalization of the previous argu- 
ments, this tensor provides in general the maximum lo- 
calizability in any given direction. An equivalent expres- 
sion for the localization tensor is: 



^Jdvj dr> (r-r>) a (r-r>)p\P(r,r>)f, 



(A7) 



which has the meaning of the second moment of the 
(squared) density matrix in the coordinate r — r'. 

At this point, we may think of a crystalline solid as 
of a very large "molecule", or a cluster, and take the 
thermodynamic limit. Since bulk properties must be 
independent of the choice of boundary conditions (either 
BvK or "free"), the density matrix and the localization 
tensor must be the same as the one previously found in 
this work. And indeed, a glance to Eq. (|l^) sho ws t hat 
it coincides with the thermodynamic limit of Eq. (A7) in 
the insulating case. As for the metallic case, our previous 
findings bear an important message concerning Boys 
localization. For a cluster of finite size, no matter how 
large, one can doubtless build localized Boys orbitals. 
But our results prove that, in the large N limit, the 
quadratic spread of these Boys orbitals diverges whenever 
the cluster is metallic. 
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FIG. 1. Convergence of the squared localization length 
with the size of the sampling grid, for the case of GaAs. We 
compare our discretized formula with the one used by MV, 
using a genuinely cubic grid: the label M' means M' xM' xM' 
within MV notations. 
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FIG. 2. Squared localization length vs. the inverse direct 
gap (theoretical and experimental), for several elemental and 
binary semiconductors. The inequality of Eq. (^) is strongly 
verified. The points corresponding to Si and Ge with the 
theoretical gaps are out of scale. 
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FIG. 3. Trends of the squared localization length vs. the 
inverse experimental minimum gap. The lines connect the 
isoelectronic series MgS-AlP-Si and ZnSe-GaAs-Ge, and the 
isovalent one C-Si-Ge. 
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FIG. 5. Hermaphrodite orbitals for GaAs. The quantity 
displayed is TO O c(a0, defined as the yz average of the square 



modulus of the orbital to, 



j , for S2 = S3 = 0, and for the 



four j values localizing within the same cell. 
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FIG. 4. Hermaphrodite orbitals for Si. The quantity 
displayed is ni oc (x), defined as the yz average of the square 



modulus of the orbital wa 



for S2 = S3 = 0, and for the 



four j values localizing within the same cell. 
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FIG. 6. Exponential decay of two hermaphrodite orbitals 
for Ge. Thin lines correspond to the logarithm of two different 
ni oc (x), defined as the yz average of the square modulus of 
the orbital Wj,s 2 ,s 3 , for the same j and two different points of 
the two-dimensional mesh(s2,S3). The one with the slowest 
decay corresponds to S2 = S3 = 0. Thick lines are the 
macroscopic average (see text) of \nn\ oc (x). The lower panel 
is a magnification of the region indicated in the upper panel 
in order to better appreciate the linear behavior of lnni oc (a;). 
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FIG. 7. Exponential decay length, averaged over 

the two-dimensional mesh (32,53), vs. our localization 
length (square root of the second cumulant moment). The 
straight-line segments are only a guide to the eye linking com- 
pounds of the same isoelectronic series. The vertical bars are 
an estimate of the accuracy of the interpolation scheme used 
to extract the b value from the asymptotic macroscopic aver- 
age of n\ oc {x). 
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